Phase diagrams of bone remodeling using a 3D stochastic cellular automaton

We propose a 3D stochastic cellular automaton model, governed by evolutionary game theory, to simulate bone remodeling dynamics. The model includes four voxel states: Formation, Quiescence, Resorption, and Environment. We simulate the Resorption and Formation processes on separate time scales to explore the parameter space and derive a phase diagram that illustrates the sensitivity of these processes to parameter changes. Combining these results, we simulate a full bone remodeling cycle. Furthermore, we show the importance of modeling small neighborhoods for studying local bone microenvironment controls. This model can guide experimental design and, in combination with other models, it could assist to further explore external impacts on bone remodeling. Consequently, this model contributes to an improved understanding of complex dynamics in bone remodeling dynamics and exploring alterations due to disease or drug treatment.


Introduction
Bone remodeling is a very complex process necessary to ensure a healthy bone structure.This process can be disturbed by a number of causes, for example hormonal imbalance, extreme load conditions or cancer-related bone disease.To understand this complex interplay, bone remodeling has been extensively studied over the last decades, including in silico models on macro-and micro-level.Early studies were focused on adaptation of the bone structure and density based on the prevailing loading conditions.Whereas first models were still run in 2D [1][2][3][4][5][6][7], later models would move to a 3D domain [8][9][10][11][12][13][14] partially include experimental data from micro-computed tomography scans [15][16][17][18][19].
These models have greatly advanced our understanding of bone architecture.However, they all focus on the (dis-)appearance of mineralized bone as a response to a mechanical stimulus.Moreover, they contain a very simple representation of the many cellular processes underlying the coupling of bone resorption and formation.Whereas the previous models have been very successful in predicting the change of bone macrostructure, there was (and still is) a need for micro-level models of multicellular interaction.Early attempts of non-spatial models were already made in the 2000s [20,21] and continued during the following decades [22][23][24][25].At the same time spatial multicellular models were introduced [26][27][28][29], allowing to illustrate the importance of local signaling mechanisms compared to global ones.More recently, some of these models have been coupled with mechanical simulations to investigate the interplay of bone resorption and formation [30][31][32].
The results of those models are promising, but most of them (except for [20,29]) are limited to specific signaling pathways.Indeed, there are still open questions and ongoing discussions about the exact interplay and origin of signaling cytokines [33][34][35][36].This discussion becomes especially relevant when bone remodeling is coupled with diseases such as cancer [37].From the perspective of in silico modeling, there is a delicate balance between adding another biochemical or mechanical pathway to better approximate the reality, and creating a model that is too complex to analyze and run, especially when complex differential equations are involved.We believe that a stochastic cellular automaton model that operates on the principles of evolutionary game theory (EGT), as introduced by Ryser et al. [29], can help to solve this contradiction.
In game theory, in its simplest formulation, each player follows one out of a set of different competing strategies.The outcome of playing one chosen strategy depends on which strategy the opponent is playing and is integrated into a matrix called payoff matrix.In EGT the payoff matrix determines the fitness of an individual based on its strategy and the relative frequency of all other strategies present in the population at a given time.The fitness eventually determines if that individual propagates and, consequently, if the relative frequency of that strategy increases or not [38].A general question in the context of EGT relates to which conditions allow different strategies to coexist.It has been found that coexistence depends on the neighboring structure of the game: spatially explicit models, in which individuals interact only with their neighbors, deliver different and broader coexisting conditions than mean-field models [29,38].
Cellular automaton models are ideal to recapitulate local dynamics of the microenvironment since they offer a straightforward definition of neighborhood.In this work, we implement a 3D stochastic cellular automaton, where voxels interact only with their nearest neighbors in a volume of interest representing bone tissue.At each time point, each voxel can take one of four different states representing the different phases of bone remodeling: Resorption, Formation, Quiescence and Environment.To create a compact representation of the frequency-dependent interaction between the local phases of bone remodeling, we consider a voxel as an individual and its state as a strategy in an evolutionary game.This representation encodes knowledge about the mutual impact the main actors of bone remodeling (osteoclasts, osteoblasts, osteocytes and environmental impacts) have on each other.Each payoff parameter in the model therefore is indirectly connected to the biological processes.Furthermore, the progression of the voxel states is explicitly related to the frequency of all strategies of a given neighborhood.In summary, EGT takes into account how bone remodeling depends both on the frequency of all strategies and on their interaction rules.
We show that a stochastic cellular automaton model governed by EGT rules is capable to simulate interactions of the microenvironment of bone remodeling by using payoff parameters instead of several differential equations.Screening single specific payoff parameters, we investigate parts of the parameter space and highlight phase transitions of the model dynamic.Associated with that, we illustrate why modeling a small neighborhood is essential for investigating local mechanisms.Assigning different payoff parameters to model physiological and pathological behavior of bone, we can draw first conclusions on the necessary interactions of Resorption, Formation, Quiescence and Environment in healthy bone.By doing so, we can identify key parameter combinations required for physiological bone remodeling, as well as the sensitivity of the system transitioning to pathological remodeling.This model opens the window to explore the complex interplay of factors essential for healthy bone remodeling and potential alterations due to disease (through local or systemic signaling, via soluble factors or hormones) or drug treatment.

Results
We created a spatial, stochastic cellular automaton model for the simulation of bone remodeling, which uses the concepts of EGT for the update criteria of the simulation domain.The four states of the cellular automaton (Resorption, Formation, Quiescence and Environment) represent sites of activity of the actors of bone remodeling (osteoclasts, osteoblasts, osteocytes and environmental impacts) (Fig 1A and 1B).Per iteration, one voxel proliferates (spreads) to one of its neighboring voxels.The initial configuration can expand during the simulation, when the initial boundaries are exceeded by a state different from Environment.Furthermore, we restrict some interactions to prevent non-physical behavior (see Materials and methods).The core of the model are the payoff parameters, which encode the impact that one voxel state has on the spreading of another one and therefore represent biological interactions (Fig 1C).The payoff parameters range from −1 to 1. Positive values stand for a stimulating impact, negative values for an inhibitory impact.For example, a Formation-voxel turns into a Quiescence-voxel only if there is a Quiescence-voxel in its immediate neighborhood.Consequently, a Quiescencevoxel needs to expand into the Formation-voxel.In this context, the Formation-voxel needs to have a positive influence on the expansion rate of the Quiescence-voxel (i.e.g QF has a positive value).For a more compact display, the payoff parameters are depicted in the payoff matrix (Fig 1D

Single patch simulations
To reduce the complexity resulting from multiple interactions, we first examine each sub-process of bone remodeling (Resorption and Formation) on its own.We set up simulations with either a Resorption-or a Formation-patch of 3x3 voxels positioned on the surface of a cube made of 28x28x28 Quiescence-voxels.In each case, the number of relevant payoff parameters reduces down to three (Fig 2).To address the stochastic nature of the model, each parameter combination runs for 20 rounds with a different random seed per round.
Resorption.When simulating the resorption process of bone remodeling (Fig 1A ), the main underlying interactions are between osteoclasts and the quiescent bone.During bone remodeling, resorption initiates (seeding of Resorption-patch), progresses for a limited amount of space and time and then stops/disappears.The main payoff parameters influencing the progression of the Resorption-state are g RR and g RQ .In addition, the speed of demineralization, that is, Environment-voxels spreading to Resorption-voxels, is denoted by g ER (Fig 2A).
We vary the parameters g RR and g RQ between −1 (maximal inhibition) and 1 (maximal stimulation) and set g ER to 0.5.We group the results in a phase diagram containing three different categories: not active Resorption (Fig 3A, gray), bounded spreading of Resorption (Fig 3A, yellow to orange) and not bounded (Fig 3A , lilac).In the not active case, the Resorption-voxels are taken over by Environment-voxels, before they can spread.In the case of bounded spreading, the Resorption-patch starts to spread and stops within the simulation time.In the case of not bounded spreading, the Resorption-patch starts to spread and does not stop within the simulation time (Fig 3C, S1 Video).
To account for stochastics, we run 20 independent simulations for each of the relevant parameter combinations.Some of the parameter combinations are either always, i.e. 20 out of 20, not active (gray) or always lead to not bounded spreading (lilac).Between those two phases, we find parameter combinations that lead to bounded Resorption in either some (< 100%) or all (100%) of the 20 rounds (Fig 3A).Whenever there is a combination of two negative parameters (bottom left corner), the result is always not active, since the impact of both parameters (g RR and g RQ ) on the Resorption patch is inhibitory.Whenever both parameters have high positive value (meaning a strong stimulation), the spreading is not bounded anymore (upper right corner).The bounded spreading phase is located between those two extremes and displays a yellow-orange gradient indicating the percentage of bounded spreading runs out of the 20 simulations.This result remains qualitatively the same, when changing the demineralization speed g ER (S1A Biologically, the case of bounded spreading is the most relevant, since this is the only dynamic that can lead to an effective but limited remodeling process.Therefore, the phase within the parameter domain leading in 100% of the simulations to bounded Resorption is the most applicable.This phase is mainly situated in the upper left corner of the phase diagram, which represents positive g RR values and negative g RQ values.This can be interpreted as a stimulating impact of osteoclast signaling on osteoclasts (g RR > 0) and an inhibitory impact of osteocyte signaling on osteoclasts (g RQ < 0).The diagonal alignment of this phase suggests that the higher the stimulation from the osteoclasts, the higher the inhibition required from the osteocytes to keep Resorption bounded.Furthermore, that phase is more stable when the absolute value of both parameters is high.Focusing on only those parameter combinations leading in 100% of the simulations to bounded Resorption, we evaluate the median of net decrease in quiescent volume (Fig 3B).Here we find, that on the border to the not bounded spreading phase, more voxels are resorbed than on the border to the not active phase.
Formation.When simulating the formation process of bone remodeling (Fig 1A ), the main underlying interactions are between osteoblasts and the quiescent bone.During bone remodeling, formation initiates (seeding of Formation-patch), progresses for a limited amount of space and time and then stops/disappears.The main payoff parameters influencing the progression of the Formation-state are parameters g FF and g FQ .In addition, the speed of mineralization, that is, Quiescence-voxels spreading to Formation-voxels, is denoted by g QF (Fig 2B).
We vary the parameters g FF and g FQ between −1 and 1 and find patterns similar to the Resorption-simulations.Again, we group the results in a phase diagram containing three differ-  Again, we run 20 independent simulations for each of the relevant parameter combinations.Some of the parameter combinations are either always, i.e. 20 out of 20, not active (gray) or always lead to not bounded spreading (green).Between those two phases, we find parameter combinations, which lead to bounded Formation in either some (<100%) or all (100%) of the 20 rounds ( The bounded Formation (100% of the simulations lead to bounded spreading) and the parameter combinations leading to it are the most relevant.They are mainly situated in the bottom right corner of the phase diagram, which represents negative g FF values and positive g FQ values.This can be interpreted as an inhibitory impact of osteoblast signaling on osteoblasts (g FF < 0) and a stimulating impact of osteocyte signaling on osteoblasts (g FQ > 0).In contrast to Resorption, the bounded Formation is not clearly aligned diagonal, but rather spans through the lower left triangle of the fourth quadrant.It is more stable for a high inhibitory impact of osteoblasts on osteoblast, and gets less stable as this inhibitory impact declines.This is further discussed in the Discussion section.Focusing on only those parameter combinations that lead to 100% bounded spreading, we evaluate the median of the net increase in quiescent volume (Fig 4B).Again, we find that on the border to the not bounded spreading phase, more voxels are formed than on the border to the not active phase.

Effect of the size of the neighborhood
The choice of the cellular automaton neighborhood has a strong influence on the overall simulation dynamics.The neighborhood in the single patch simulations is initially set to a von-Neumann neighborhood (radius = 1), meaning that only voxels sharing a surface with the voxel in question are considered neighbors (Fig 5A).For a 3D grid of cubic cells, this involves six voxels.To understand the role of the size of the microenvironment on the model dynamics, the single patch simulations are repeated with the next larger neighborhoods also typical for cellular automaton models.For this, we choose a Moore neighborhood, meaning that all voxels sharing a surface or a corner with the voxel in question are considered neighbors.Two sizes of Moore neighborhood are considered: radius = 1 with 26 voxels on a cubic grid and radius = 2 with 124 (Fig 5B and 5C).Finally, a mean field approximation, where all voxels in the simulation are considered neighbors, represents the largest neighborhood possible (Fig 5D ).
When increasing the size of the neighborhood only slightly beyond the nearest neighbors (Moore r = 1, Fig 5B), the phase of bounded spreading in both cases (Resorption and

The remodeling cycle: Combining Resorption and Formation
To simulate one cycle of bone remodeling, we now use the results from the single patch simulations, focusing on the parameter sets that lead in 100% of rounds to bounded spreading.For this, we choose different combinations of speed of demineralization (g ER ) and speed of mineralization (g QF ).Then we combine all parameter sets for Formation (g FQ , g FF ) and Resorption (g RQ , g RR ), which lead in 20/20 runs to a bounded spreading phase for the respective speed.Both processes (Resorption and Formation) are seeded one after the other, thereby assuming that the two processes are separated in time.Since bone remodeling is the process of bone maintenance and homeostasis, we monitor the net change of volume at the end of the remodeling cycle.
To imitate the pattern of the bone remodeling cycle from Fig 1A, we re-fill some of the Environment-voxels left behind by Resorption with Formation-voxels.An example of this is shown in Fig 7A, where we combine two parameter pairs leading to bounded spreading in Formation (g FF = −0.85,g FQ = 0.4) after Resorption (g RR = 0.9, g RQ = −0.6)for g QF = 0.3 and g ER = 0.8 to simulate one remodeling cycle (see also S3 Video.After simulating all combinations of parameter pairs (g RR , g RQ ) and (g FF , g FQ ) for the three combinations of (de-)mineralization speed (g ER < g QF , g ER = g QF and g ER > g QF ), we then quantify the total volume change at the end of the remodeling cycle.
In contrast to the single patch simulations, some parameter combinations lead to unbounded Formation in individual runs.This affects less than 2.5% of total parameter combinations for each of the three cases (for details, see S1 Table ).We excluded all parameter combinations with unbounded Formation runs for the evaluation in Fig 7B .The distribution of the net increase in volume, however, differs only slightly compared to the evaluation of all In all three combinations of (de-)mineralization speed, the median of the total volume change is slightly larger than zero (Fig 7B).This implies a net volume growth at the end of a remodeling cycle.This tendency increases with the ratio of speed of demineralization to mineralization, g ER g QF .

Discussion
With our cellular automaton model, we can explore the complex interactions necessary for a well-balanced bone remodeling cycle.Although the different states of bone remodeling (Resorption, Formation, and Quiescence) are not based on an evolutionary process in the classical sense, the prevalence of one state is directly related to the frequency of the other two and itself.Therefore, EGT is well suited to theoretically explore the dynamics in the bone microenvironment that lead to perturbations in homeostasis, as it takes into account the dependence on the frequency of all strategies as well as their interactions.Using the concepts of EGT enabled us to abstract various signaling pathways with just two payoff parameters per interaction, without making model assumptions about complex in vivo signaling cascades.Implementing EGT as a cellular automaton allowed us to highlight the importance of local dynamics in the bone microenvironment with respect to the spread of the Resorption and Formation state.In view of future applications of the model (for example to study interactions between bone remodeling and cancer cells), EGT makes it easy to expand the model without major changes of the model philosophy and its implementation.For each newly introduced strategy, only the payoff matrix would need to be expanded with additional entries.However, simplicity also has its limitations: Parameterizing EGT models with experimental data is difficult, as it is for computational models of complex biological behavior in general.Although the payoff parameters cannot be traced back to specific signaling pathways, much can be learned about the interactions of the bone remodeling actors and how their frequency,

Resorption
With the phase diagram of the single patch Resorption simulations we showed how changing either one of the studied parameters (g RQ or g RR ) can change the whole dynamic of the resorption process: from inactive Resorption over bounded to unbounded Resorption (Fig 3A).We considered the phase of bounded Resorption as physiological behavior, since Resorption progresses first, but stops at a certain point.We further observed variation in the amount of resorbed bone within the bounded spreading phase (Fig 3B).This, again, was a result of varying one or both of the studied parameters.
When we investigate the parameter combinations leading to bounded Resorption, we can put this in context with current discussions on bone remodeling signaling.Biologically, Resorption is mainly stimulated through macrophage colony stimulating factor (MCSF) and receptor activator of nuclear factor κB (RANK) originated from cells of the osteoblast lineage [35].Those, however, should not be confused with active osteoblasts involved in bone formation, but rather include osteocytes within the bone matrix, bone lining cells or osteoblast progenitors [35,36].Consequently, this stimulation is represented by the parameter g RQ .Furthermore, parameter g RR represents autocrine regulation of osteoclast formation by cytokines interleukin-1α (IL-1α) and tumor necrosis factor α (TNF-α), as well as osteoclast progenitors recruitment regulated by sphingosine-1-phosphate (S1P) derived from osteoclasts [33,36].
When both parameters (g RQ and g RR ) are positive (therefore stimulating) the model leads to unbounded spreading of Resorption (Fig 3A).For bounded spreading, one of the two stimuli needs to be negative.Moreover, the bounded Resorption phase is more stable for a strong g RR stimulus than for a strong g RQ stimulus.The fact that all of the 100% bounded spreading runs are located at positive g RR and negative g RQ stimulus, suggests that living osteocytes negatively regulate osteoclasts [35].According to this, the recruitment of osteoclasts is caused by osteocyte death (due to microcracks, low oestrogen levels, etc.) and the consequently weakened (meaning less negative g RQ ) regulation of osteoclasts.Possible spikes of stimulation by osteocyte death to initiate Resorption (recruit osteoclasts [39,40]) are not yet captured by the current model.Here, we simulate the initiation by seeding a patch of Resorption voxels in the initial configuration.

Formation
With the phase diagram we obtained from the parameter analysis of the Formation single patch simulations, we showed how changing either one of the studied parameters (g FQ or g FF ) can change the whole dynamic of the formation process similar to Resorption.We found the same three phases: inactive, bounded and unbounded Formation (Fig 4A).Again, we considered the phase of bounded Formation as physiological behavior and observed variations in the amount of formed bone within the bounded spreading phase (Fig 4B).The differences between the shape of the bounded spreading of Resorption (Fig 3A ) and the bounded growth of Formation (Fig 4A ) arose from the fact that Resorption spread into the Quiescence-voxels, which had an impact on the seeded Resorption/Formation-voxels, whereas Formation spread into the Environment-voxels, which had no impact (yet).Therefore, Formation was more tolerant about a stimulating impact from Quiescence (g FQ ) than Resorption (g RQ ).
When we investigate the parameter combinations leading to bounded Formation, we can link this with current discussions on bone remodeling signaling.For formation, osteoblasts are stimulated and/or inhibited through bone matrix either by osteocytes (stimulation through messengers, inhibition through sclerostin) [35] or osteocalcin (inhibition) [34] and biglycan (stimulation) [34].This is represented by the parameter g FQ .The parameter g FF describes the impact of Formation on Formation and includes but is not limited to the release of osteoblastderived insulin-like growth factor 1 (IGF-1), which leads to osteoblast differentiation [36].Additional stimulation through growth factors released during demineralization of the bone matrix is not represented for now.This concerns the impact of Environment on Formation (g FE ), but would be limited to Environment voxels that were recently generated due to Resorption.Technically, this would be possible by tracking the age of each voxel.For now, however, the purpose is to keep the model as simple as possible, which includes keeping constant values for the payoff parameters in space and time and assuming the impact of Environment (g iE ) to be zero.
Similar to the results for the Resorption patch, the spreading of Formation is unbounded when both parameters (g FF and g FQ ) are stimulating (Fig 4A).Most of the 100% bounded spreading runs are located at positive g FQ and negative g FF stimulus, which suggests that stimulation through the bone matrix is the most important mechanism [35].Cases of bounded spreading can also be found for positive g FF and negative g FQ , although never for 20 runs.This suggests that, for bounded Formation, positive stimulation through autocrine factors can to be compensated by negative regulation through the bone matrix, but this regulation would not be stable.To our knowledge, autocrine inhibition (therefore g FF < 0) of osteoblasts is not discussed in the literature, yet.

Effect of the size of neighborhood
Comparison of different neighborhood sizes demonstrated that spatial modeling with nearest neighbors can reveal local control mechanisms of the microenvironment that are otherwise lost in the mean field approximation (Fig 6).This agreed with previous findings of [29], who showed for certain parameters spaces that an unstable EGT dynamic can stabilize when switching from a non-spatial to a spatial model.
The size of the bounded Resorption phase became narrower with increasing size of the neighborhood (Fig 6).The same applied for Formation (Fig 6).Furthermore, comparing neighborhood sizes showed, that differences between the shape of bounded Resorption (Fig 3A ) and bounded Formation (Fig 4A ) in a 6-voxel-neighborhood reduced progressively with an increasing neighborhood.This confirmed, that those differences arose from different voxeltypes, and therefore directions, in which Resorption and Formation spread: the former into the Quiescence-voxels, the latter into the Environment-voxels.This showed how regulation mechanisms that work locally within a small microenvironment can be overlooked, when modeling with too large a neighborhood or a mean field approximation.

The remodeling cycle: Combining Resorption and Formation
The combination of bounded-spreading Resorption and Formation gave promising first results of simulating a whole bone remodeling cycle using parameters from the bounded spreading phase.Quantitatively, the results were still tilted towards bone growth instead of preserving the bone volume (maintenance).This suggested that if the processes of Formation and Resorption would be repeated over and over again (simulating a constant remodeling behavior), the simulated domain would most likely increase its volume over time.This however does not mirror the maintaining nature of the bone remodeling process.
The fact that some parameter combinations for Formation always led to bounded spreading in the single patch simulations, but not when combining Resorption and Formation, emphasized again the importance of the local environment.When simulating the remodeling cycle, Resorption acted as anticipated from the results of the single patch simulations, because it was seeded in the same way (3x3 voxel-patch on the flat surface).Formation, however, was seeded differently: Instead of seeding on a flat surface, we based the spatial distribution of the Formation voxels on the Environment-pattern that Resorption left behind.This caused in some cases differences in behavior of Formation compared to the single patch simulations.

Perspective
In Figs 3B and 4B we show that different parameter combinations within the phase of bounded spreading lead to different amounts of net decrease and increase in quiescent volume, respectively.To ensure maintenance of the volume, the parameter pairs for Resorption (g RR and g RQ ) and Formation (g FF and g FQ ) might need to be adjusted to each other.This could lead to a subphase within bounded spreading, that results in not just bounded but in physiological bone remodeling.Alternatively, some of these parameters (e.g.g FF and g FQ ) could be dependent on the values of others (e.g.g RR and g RQ ).
Criteria for distinguishing physiological bone remodeling from pathological bone remodeling may be deduced from the experimental investigations of Young and colleagues [41,42].The spatially explicit nature of our model enables the integration of bone geometry, as characterized in the experimental data set, to serve as the initial configuration for the model.This should allow for a direct juxtaposition of the bone remodeling patterns observed experimentally with those forecasted by the model.As a first approach, one could use the quantification of resorbed and formed bone volume over time and the total change in quiescent bone volume from the experimental data to identify payoff parameter ranges that are biologically meaningful.This parameterization could be based on the results presented in this work.Additionally, comparisons with experimental data from different species would be possible without changing the model, but taking into account potential species-specific differences of the identified payoff parameter ranges.Furthermore, interactions that are found to be of particular interest could be additionally modeled by replacing the constant value of the payoff parameter with equations of a submodel that determine their time or spatial dependence.
Nevertheless, to effectively implement this approach, a constellation of technical and theoretical challenges must be addressed: the temporal dimension of the simulations must be harmonized with the biological tempo of bone remodeling by calibrating the rates within the payoff matrix; the biological variation between individual mice requires an evaluation of the margin of error within which the model reliably mirrors the experimental findings; and, while the experimental data captures alterations in quiescent bone, the model simulates voxel dynamics over time.We remain optimistic that forthcoming research endeavors will resolve these complexities.
So far, the impact of the Environment state is kept neutral (without any impact) and consequently the parameters g iE are set to zero.However, systemic regulation by hormones and local regulation via soluble factors modulate bone physiology.Alterations of these as well as the effect of drugs could be integrated and modelled through the Environment-state.Concretely, this means to investigate the influence of the parameters g iE on the previously discussed phase diagrams.
For now, the number of Formation voxels being seeded after resorption is the same as the previously seeded Resorption voxels.The legitimization of that assumption can be further investigated by varying the number of seeded Formation voxels and possibly coupling that number to the size of the cavity left behind by Resorption (Fig 7A, upper right corner).
Since both, Resorption and Formation patch, are seeded on different timescales, the impact that Resorption and Formation have on each other (g RF and g FR ) is not captured in the current setup.This assumption is supported by Sims and Martin [36] who showed that coupling between active osteoclasts and osteoblasts via direct contact is unlikely.Findings of [43], however, introduce a mixed reversal-resorption phase, where at the end of a resorption phase the density of osteoprogenitor cells grows until a certain threshold is reached and formation starts.This dynamic could be modeled by seeding the Formation patch while Resorption voxels are still present, thereby simulating the switch from reversal-resorption-to formation-phase by Formation voxels out competing remaining Resorption voxels.
In the proposed model, the Resorption and Formation patch are seeded manually.In the future, the model could be coupled with a finite element analysis similar to [12].The resulting strains would inform the seeding of the remodeling patches.According to this, the Resorption patches would be seeded randomly, assuming regular appearing microcracks and other damages in the bone.The Formation patch, however, would be seeded in areas of high strains, which can either be caused by increased load conditions or a resulting resorption cavity.Furthermore, the magnitude of strain could also serve as a criterion for the number of Formation voxels being seeded.

Conclusion
With the proposed model, we laid the foundations for simulating spatial bone remodeling with a stochastic cellular automaton using evolutionary game theory.Such modeling approach allowed us to generate phase diagrams of bone resorption and formation to qualitatively understand the impact and sensitivity of stimulatory/inhibitory effects of certain biological processes.Furthermore, comparing the effect of different neighborhood sizes demonstrated that spatial modeling with nearest neighbors can reveal local control mechanisms of the microenvironment that are lost in the mean field approximation.
In the future, the proposed model could inform in vitro/in vivo experimental design and/or be combined with other in silico models investigating external impacts on bone remodeling (finite element analysis of mechanical loading, cancer models or drug/treatment protocols).To sum up, this model opens the window to explore the complex interplay of factors essential for healthy bone remodeling and potential alterations due to disease (through local or systemic signaling, via soluble factors or hormones) or drug treatment.

Model
A spatial four strategy model resembling a stochastic cellular automaton, whose update rule is based on criteria borrowed from EGT is set up in Matlab (version 2023b).The simulation domain is approximated by a cube and discretized into n same-sized-voxels per axis.Each voxel gets assigned a voxel state.Based on the process of bone remodeling, the voxel states can vary between Environment (= 1), Resorption (= 2), Formation (= 3) and Quiescence (= 4) (Fig 1A and 1B).The current state of a voxel with position x is described by ξ(x) 2 {1, 2, 3, 4}.We adopt the notation of Ryser et al. [29], where the payoff parameter g ij describes the impact each state j has on the spreading rate of the state i (Fig 1C).The set of voxels previously defined as neighbors of voxel x is denoted by y * x.The neighborhood is set to include only the nearest neighbors, meaning those which share a surface with the voxel of interest (Fig 5A ).
All interactions described with payoff parameters g ij are coded in a payoff matrix (Fig 1D), which contains sixteen entries for a system with four strategies.They are decisive for the dynamics of the whole system.The payoff parameters range between -1 (maximally inhibitory) and 1 (maximally stimulating) to capture their relative proportion to each other.The payoff parameters are constant in space and time and are defined at the beginning of the simulation.
The dynamics of the simulated system are calculated by applying a Gillespie algorithm, which uses the expansion rates of the voxels to determine their probability to spread.The final expansion rate for each voxel per iteration is calculated by summing up the impact from all neighboring voxels y: After the spreading rates C(x) for all voxels have been calculated, each voxel x gets assigned a random time exponentially distributed with mean time C(x) −1 (with C(x) −1 = 1 if C(x) = 0).The voxel that gets assigned the shortest time spreads into one of the k voxels within the neighborhood, chosen uniformly at random.
Before spreading, the neighborhood states of this voxel are checked for non-physical changes (Table 1), which are then excluded (e.g.Formation spreading to Quiescence or Quiescence spreading to Resorption).The decision whether the spreading from Resorption to Formation and vice versa are to be excluded or not is up for discussion and needs further investigation.For this work, we model Resorption and Formation on separated time scales.Consequently, they are not simultaneously present in the model domain and our results are independent of that decision.
As boundary condition, we assume Environment-voxels outside of the simulated domain, which partially constitute the neighborhood of the boundary voxels.Each time the initial boundaries are exceeded by a state different from Environment, the domain expands by one layer of voxels on each side of the cube.This enables a growing and shrinking simulation domain.

Single patch simulations
Due to the high number of payoff parameters, it is difficult to investigate all of them at the same time.We assume the impact of the Environment-voxels to be zero for this work (all parameters g iE = 0, yet in the future this could include the impact of disease and drugs.Furthermore, we do not expect any interaction connected to spreading between Quiescence and Environment (g EQ = 0, g QE = 0).
To reduce the parameters further, we set up simulations of only Resorption or only Formation, i.e. single patch simulations.For this, we seed a 3x3 voxel Resorption-or Formation-patch on the surface of a 28x28x28 voxel Quiescence-cube, which is surrounded by one layer of Environment-voxels.The whole domain has a size of 30x30x30 voxels.
In each case, this reduces the number of relevant parameters down to three.First, when simulating the single Resorption-patch, no Formation is present or can emerge and therefore all parameters g iF and g Fi are zero.Also, without Formation present, Quiescence will not be able to spread beyond its initial voxels, so we can assume parameter g QR as zero as well.Second, the argumentation line for the single Formation-patch simulation is similar.All parameters g iR and g Ri are zero.Also, without Resorption present, Environment will not be able to spread beyond its initial voxels, so we can assume parameter g EF as zero as well.
In the end, the single Resorption-patch simulation is modeled with the payoff parameters g RR , g RQ and g ER and the single Formation-patch simulation is modeled with the payoff parameters g FF , g FQ and g QF .
The parameter g ER represents the rate of demineralization, does not have a direct impact on the spreading of Resorption and is assigned a fixed value of 0.5 (effect of values 0.3 and 0.8 in S1A Fig) .Analogously, the parameter g QF represents the rate of mineralization, does not have a direct impact on the spreading of Formation and is assigned a fixed value of 0.5 (effect of values 0.3 and 0.8 in S1B Fig) .Finally, we vary g RR and g RQ between -1 and 1 in the single Resorption-patch simulation, while we vary g FF and g FQ between -1 and 1 in the single Formationpatch simulation.
Considering the stochastics of the model, we run 20 simulations (with 20 different random seeds) for each parameter combination.Each simulation runs for 1000 iterations.We expect some occasional changes (unbounded runs becoming bounded) for a significantly larger amount of iterations.However, those occasional cases do not change the main resulting phase diagram and therefore have no influence on the conclusion.
The results are categorized in the following way: Simulations, where the patch has disappeared without initial spreading, are considered not active.When the patch has spread, but still disappeared at the end of the simulation, it is considered as bounded spreading.Simulations, where the patch has not disappeared by the end of the 1000 iterations, are considered as unbounded spreading.Depending on the parameter combination, all 20 simulation runs can either fall in the same category or can fluctuate between different ones.

Effect of the size of neighborhood
To investigate the impact of the size of the neighborhood on the model dynamics, we repeat the single patch simulations (20 runs a ´1000 iterations with g ER = 0.5 and g QF = 0.5) for four different neighborhood sizes: von-Neumann neighborhood with radius 1 (6 neighbors in 3D), Moore neighborhood with r = 1 (26 neighbors in 3D) and r = 2 (124 neighbors in 3D)and mean field approximation, where all voxels in the simulation domain are considered neighbors (Fig 5).The neighborhood specifies a) which neighboring voxels are considered for the calculation of the spreading rates, and b) where a voxel can theoretically spread to.
In the case of the mean field approximation, each simulation is run for one million iterations.Furthermore, runs are only considered bounded, if less than 10% of the Quiescence-volume have been resorbed or formed.

The remodeling cycle: Combining Resorption and Formation
To simulate one whole bone remodeling cycle, one single Resorption-patch simulation is followed by one single Formation-patch simulation.We only use parameter combinations (g RR and g RQ ; g FF and g FQ ), which lead to a bounded spreading in all runs of the single patch simulations.To investigate the relation of g ER and g QF , we simulate three different combinations: g ER = 0.3 and g QF = 0.8, g ER = 0.5 and g QF = 0.5, and g ER = 0.8 and g QF = 0.3.In all three cases, we combine all parameters sets that lead to bounded Resorption (g RR and g RQ ) and bounded for Formation (g FF and g FQ ).Again, we run 20 simulations (with 20 different random seeds) for each parameter combination.
The two processes of Resorption and Formation are combined on separated timescales.First, a Resorption patch of 3x3 voxel is set on the surface of a 28x28x28 quiescent cube (identical to the single patch simulation, Fig 2A).Since we only take the parameter combinations for bounded spreading, we can be sure that after 1000 iterations all Resorption voxels will have disappeared from the model domain.Next, we set nine Formation voxels in the deepest points of the cavity that Resorption has left behind (Fig 7A, lower right).Due to the changed Formation setup compared to the single patch simulation, we can not expect that all Formation-voxels will have disappeared after 1000 iterations.Therefore, the simulation runs again until all Formation voxels have disappeared or the model domain exceeds a side length of 46 voxels.In case, the simulation run stops because of the latter, we exclude the respective parameter combinations (all 20 runs) from the evaluation (for details, see S1 Table ).
Fig), the size of the Resorption patch (S2A Fig) or the position of the Resorption-patch (S3A Fig).

Fig 1 .
Fig 1. Relation of the cellular automaton model and the biology of bone remodeling.A: During bone remodeling, old bone gets resorbed by osteoclasts, followed by new bone being formed by osteoblasts for the purpose of maintenance.Mechanosensitive osteocytes located in the quiescent bone mediate that process.B: The four states of the cellular automaton represent the activity of the actors of bone remodeling.C, D: The payoff parameters represent the impact that one voxel state has on the spreading of another one.They summarize all biological interactions taking place between two states and can be interpreted as stimulating (positive value) or inhibiting (negative value).https://doi.org/10.1371/journal.pone.0304694.g001 ent categories: not active Formation (Fig 4A, gray), bounded spreading of Formation (Fig 4A, yellow to orange) and not bounded (Fig 4A, green).In the not active case, the Formation-voxels are taken over by Quiescence before they can spread (Fig 4C, S2 Video).

Fig 2 .
Fig 2. Setup for single patch simulations.A: Initial configuration for the single patch simulation of Resorption (lilac) on a Quiescence-cube, where the relevant payoff parameters are g ER , g RR and g RQ .B: Initial configuration for the single patch simulation of Formation (green) on a Quiescence-cube, where the relevant payoff parameters are g QF , g FF and g FQ .All parameters marked in gray are assumed to be zero.The parameters marked in white have a constant value, while the parameters marked in orange are screened from -1 to 1. https://doi.org/10.1371/journal.pone.0304694.g002 Fig 4A).Areas where Formation is never active, are situated in the bottom left corner, while areas where Formation is never bounded, are situated in the upper right corner (Fig 4A).The bounded spreading phase is located between those two extremes and displays a yellow-orange gradient indicating the percentage of bounded spreading runs out of the 20 simulations.The three different phases also remain for Formation, when changing the mineralization speed g QF (S1B Fig), the size of the Formation patch (S2B Fig) or its position (S3B Fig).

Fig 3 .
Fig 3. Results for single patch simulations: Resorption.A, C: Summary of all simulations over the combination of parameters g RR and g RQ and spreading activity of the seeded Resorption-patch.Simulation results are grouped in always not active (gray) or always not bounded spreading (lilac).Between those two phases, we find parameter combinations, which lead to bounded spreading in either some (< 100%) or all (100%) of the 20 simulation runs.B: Median of final resorbed volume for parameter combinations that lead to bounded spreading in all runs.https://doi.org/10.1371/journal.pone.0304694.g003

Fig 4 .
Fig 4. Results for single patch simulations: Formation.A, C: Summary of all simulations over the combination of parameters g FF and g FQ and spreading activity of the seeded Formation-patch.Simulation results are grouped in always not active (gray) or always not bounded spreading (green).Between those two phases we find parameter combinations, which lead to bounded spreading in either some (<100%) or all (100%) of the 20 simulation runs.B: Median of final formed volume for parameter combinations that lead to bounded spreading in all runs.https://doi.org/10.1371/journal.pone.0304694.g004